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Abstract. - We solve the Dyson equation for atoms and diatomic molecules within the GW 
approximation, in order to elucidate the effects of self-consistency on the total energies and 
ionization potentials. We find GW to produce accurate energy differences although the self- 
consistent total energies differ significantly from the exact values. Total energies obtained from 
the Luttinger-Ward functional _Elw[G] with simple, approximate Green functions as input, are 
shown to be in excellent agreement with the self-consistent results. This demonstrates that the 
Luttinger-Ward functional is a reliable method for testing the merits of different self-energy 
approximations without the need to solve the Dyson equation self-consistently. Self-consistent 
GW ionization potentials are calculated from the Extended Koopmans Theorem, and shown to 
be in good agreement with the experimental results. We also find the self-consistent ionization 
potentials to be often better than the non-self-consistent Go Wo values. We conclude that 
GW calculations should be done self-consistently in order to obtain physically meaningful and 
unambiguous energy differences. 



Introduction. - Green function methods have been used with great success to calculate a 
wide variety of properties of electronic systems, ranging from atoms and molecules to solids. 
One of the most successful and widespread methods has been the GW approximation (GWA) 
[1], which has produced excellent results for band gaps and spectral properties of solids [2, 
3], but so far has not been explored much for atoms and molecules, although it has been 
known that for atoms the core-valence interactions are described much more accurately by 
GW than Hartree-Fock (HF) [4]. Moreover, the GW calculations are rarely carried out in 
a self-consistent manner, and the effect of self-consistency is for this reason still a topic of 
considerable debate [5,6]. In this paper we present self-consistent all-electron GW (SC-GW) 
calculations for atoms and diatomic molecules. The reason for doing these calculations is 
two-fold: Firstly we want to study the importance of self-consistency within the GW scheme. 
Such calculations are usually avoided due to the rather large computational effort involved. 
It has been suggested that self-consistency will in fact worsen the spectral properties, though 
calculations on silicon and germanium crystals indicate that this is not always the case [5]. 
The second reason is that we aim to study transport through large molecules and molecular 
chains, where it is essential to account for the screening of the long range of the Coulomb 
interaction. The calculations on diatomic molecules are the first step in this direction. 

The GWA is obtained by replacing the bare Coulomb interaction v in the exchange self- 
energy with the dynamically screened interaction W, such that £ = —GW. The screened 
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interaction also depends on the Green function, and one thus needs to solve a set of coupled 
equations for G and W . One usually goes through only a single iteration of this scheme. 
With an initial Green function Go calculated from, e.g., the local density approximation 
(LDA), one calculates W and E, and subsequently obtains a new Green function from the 
Dyson equation. This scheme, known as the Go Wo approximation, has produced good results 
for a wide variety of systems [2], but suffers from a dependence on the choice of the initial 
Go- Moreover, observables like the total energy are not unambiguously defined, and can be 
calculated in several different ways. These problems can be cured by performing self-consistent 
calculations [7], since the GWA is a ^-derivable approximation (see Fig.^l. The fact that self- 
consistency removes these ambiguities does not imply that the results are necessarily closer 
to the exact values. For the electron gas it was shown that self-consistency actually worsens 
the spectral properties, while the total energy is in excellent agreement with Monte-Carlo 
results [8]. On the other hand, for a system of very localized interactions, SG-GW produced 
poor results for both total energies and spectral properties [9]. Furthermore, Delaney et 
al. [6] recently published SG-GW results for the ionization potential of the Be atom that were 
worse than those of GoWV Calculations on the Si and Ge crystals have, however, shown that 
self-consistency leads to improved band gaps [5]. 

General formulation. In this paper, we study the importance of self-consistency in 
GW for atoms and diatomic molecules. We compare the self-consistent total energies to those 
obtained from the Luttinger-Ward (LW) functional [10] which was earlier used to estimate 
the GW total energy for atoms [11] and the electron gas [12]. The LW functional E\^[G] is 
a variational energy functional in the sense that 5E\jw[G]/5G = 0, when G is a self-consistent 
solution of the Dyson equation. This variational property suggests that evaluating -Elw on 
an approximate Green function obtained from, e.g., HF or LDA calculations will give a result 
very close to the self-consistent value. This was earlier shown to be the case for the second- 
order self-energy [13], and investigating the stability of the LW functional also for the GWA 
is an important goal of this paper. The previously published LW calculations [11] indicated 
that the GW total energies are not very accurate, but the essential question is rather whether 
total energy differences are produced accurately. We have for this reason also calculated the 
binding curve of the H2 molecule and two-electron removal energies /S.E — En— 2 — ^n- 

We use the finite temperature formalism, with a temperature T (we are only considering 
the limit T — > 0) and a chemical potential fi. The Green function depends on the imaginary 
time coordinate r, in the range —j3 < r < (3 = 1/ksT, where ks is the Boltzmann constant. 
It satisfies the Dyson equation 

V 2 

- d T + — - w(r) - v H (r) + n 

= 5(T)5(x-x') + y^dr 1 J dx 1 E[G](x,x 1 ;r-r 1 ) xG(x 1 ,x';r 1 ), (1) 

where x — (r,cr) denotes the space- and spin coordinates, w(r) is the external potential, 
E[G](x, x';t) is the self-energy and t)ff(r) is the Hartree potential. The last two objects 
are functionals of the Green function, and the Dyson equation should therefore be solved 
self-consistently, together with the boundary conditions G(x, x',r — 0) = — G(x, x';r) and 
G(x, x'; 0+) - G(x, x'; 0") = -<5(x - x'). 

In the GWA (Fig. ^) the electronic self-energy is given by E = —GW using the screened 
interaction W = v + vPW, where v is the bare Coulomb interaction l/|r — r'| and P = GG 



G(x,x';r) = 
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Fig. 1 - The GW self-energy E is the functional derivative of a functional $[G]. 

is the polarizability [1]. The Green function is transformed into a r-dependent matrix by 
expanding it in a basis of molecular orbitals obtained from an initial HF calculation. These 
molecular orbitals are linear combinations of Slater functions located on the atomic centers 
[14]. The Green function, the S[G] and the W are peaked around the endpoints (r = and 
r = ±/3) [5, 13] so their representation on an even-spaced grid is inconvenient. Instead, we 
used a mesh which is dense around the end points [5]. 

Since we calculate the Green function on the imaginary time axis, it is inconvenient to 
calculate the ionization potentials by finding the poles of the Green function in frequency 
space, G(w). We have instead used the extended Koopmans theorem (EKT) [15] where the 
ionization potentials are found from the eigenvalue equation 



where Ay = — 9 r G,j(r)| T= o, the density matrix is given by pij = Gij(0~) and the matrix 
indices refer to the molecular orbital basis [13]. The eigenvalues A m are interpreted as A m = 
E^~ x — Eq + p,, i.e. the ionization potentials plus the chemical potential. The EKT is known 
to be exact for the lowest ionization energies, if the exact A and p matrices are given [16]. 
For the HF approximation, the EKT eigenvalues obviously agree with the poles of the HF 
Green function, and it is an unproven conjecture that these two methods will give the same 
value for the first ionization potential when the Green function is calculated self-consistently 
within a conserving approximation. The EKT has recently been used to calculate ionization 
potentials for atoms and molecules from a self-consistent Green function using the second 
order diagrams [13]. 

To calculate the SC-GW total energy E = T + V ne + Uo + U xc , we use the fact that the 
exchange-correlation part of the interaction energy is given by 



and the kinetic energy T, nuclear-electron attraction energy V ne and Hartree energy Uq are 
trivially obtained from the density matrix p. There are many other ways to calculate the 
total energy from a given Green function, but only for a self-consistent solution of the Dyson 
equation will these methods give the same result [7] . One alternative is to calculate the energy 
from variational functionals of the Green function. LW have shown [10] that the total energy 
can be written as 




(2) 





E h w[G] = *[G] - U - Tr{EG} - Trln[S - G^] + /iN 



(4) 
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Table I - Total energies (in Hartrees) calculated from SC-GW compared to CI values and results 
from the LW functional and Galitskii-Migdal formula evaluated on Ghf ■ 



System 




Etw [Ghf] 


-Bgm[Ghp] 


CI 


He 


-2.9278 


-2.9277 


-2.9354 


-2.9037 1 


Be 


-14.7024 


-14.7017 


-14.7405 


-14.6674 1 


Be 2 + 


-13.6885 


-13.6885 


-13.6929 


-13.6556 1 


Ne 


-129.0499 


-129.0492 


-129.0885 


-128.9376 1 


Mg 


-200.1762 


-200.1752 


-200.2924 


-200.053 1 


Mg 2 + 


-199.3457 


-199.3453 


-199.3785 


-199.2204 1 


H 2 


-1.1887 


-1.1888 


-1.1985 


-1.133 2 


LiH 


-8.0995 


-8.0997 


-8.1113 


-8.040 3 


'From Ref. [17]. 


2 From Ref. [18]. 


3 From Ref. [19]. 







where Gh is the Hartree Green function, and £ = 8&/5G. The trace indicates an integration 
over the spatial coordinates and r [11], see also Eq. ©. It is easily verified that 5Eyw/8G = 
when G is a self-consistent solution of the Dyson equation Q). Hence, if we evaluate the LW 
functional on a simple input Green function, we obtain a result close to the self-consistent 
energy, since we make an error only to second-order in the deviation from the self-consistent 
G. This means that we have a computationally cheap way of obtaining self-consistent total 
energies. The quality of the energies will ultimately be determined by the chosen self-energy 
approximation. 

Within a molecular orbital basis, the Dyson equation JQ) becomes a matrix equation. We 
introduce a reference Green function Go in order to write the equation on integral form, 

Gij(r) = SijG ,i{T) + / dn d7 2! y)Go,i(T-Ti)E tt (r 1 -7a)Gfci(7a), (5) 
Jo Jo k 

where X = £[G] — S , and E is the self-energy corresponding to G [13]. We take G and S 
to be the HF Green function and self-energy, but this choice is arbitrary. Using, e.g., LDA 
instead would not change any of the results. The inverse temperature is chosen to have a 
sufficiently large value, typically larger than 100 a.u.. The value of the chemical potential is 
somewhat arbitrary, but should be in the gap between the highest occupied and the lowest 
unoccupied orbital. We checked that the observables calculated from the resulting Green 
function did not depend on the choice of (3 and fi. The calculations on the molecules were 
done at the experimental bond lengths. 

Results. In Table [Q we show the SC-GW total energies of some atoms and small 
molecules. We have also included the -ElwIGhf] results, which are in spectacular agreement 
with the SC-GW values. This agreement is independent of the chosen basis set, and was earlier 
observed also for the second-order diagrams [13]. The third column shows the total energy 
calculated from Ghf using the Galitskii-Migdal [20] formula. In contrast to the LW results, 
these are not in good agreement with the self-consistent energies. This clearly demonstrates 
that different total energy functionals will not produce the same results when evaluated on a 
non-selfconsistent Green function (in this case, Ghf)> and it also demonstrates the importance 
of using the variational functionals for obtaining a result in agreement with the self-consistent 
values. 

As a further test of the total energy functionals, we have calculated the total energy of 
the H 2 molecule for a range of internuclear separations. Figure |21 shows the SC-GW results 



A. Stan et al: Fully SC-GW calculations for atoms and molecules 



5 



-0.8 
-0.S5 
-0.9 
7 -0.95 
I -1 
-1.05 
-1.1 
-1.15 



1 2 3 4 5 6 7 8 

R (a.u.) 

Fig. 2 - The total energy of the Kb molecule, as function of the interatomic distance, calculated within 
the second order, the self-consistent GWA, the E§w [GW] functional and CI (from Ref. [18]). For 
comparison, the HF results are also presented. 

together with the -Elw[Ghf] energy. The curves agree closely up to R « 5 a.u., and the 
deviation remains small even at R = 8. The gradual increase in the deviation is due to the 
fact that the input Ghf differs increasingly from the self-consistent Green function at large 
separations, making the variational property of -Elw l ess reliable. We also plotted benchmark 
configuration-interaction (CI) results and the binding curve obtained from the self-consistent 
Green function within the second-order self-energy approximation [13], which we were able 
to calculate up to R = 6. The second-order results are closer to the exact results than the 
GW curve around the equilibrium distance. This was to be expected, since the main feature 
of GWA is to screen the long range interactions. For atoms or small molecules it is more 
important to take both direct and exchange diagrams into account to the same order. Also 
for the atoms, the SG-GW results are not particularly close to the CI results, as seen in Table 
J] It should be noticed, however, that the shapes of the GW and the second-order curves are 
similar to each other and to the CI curve around the equilibrium bond distance. We finally 
note that, like the HF method, self-consistent GW is not a size-consistent method, i.e. the 
total energy calculated at large separations will not converge to the sum of the total energy 
of the fragments. This is not surprising, since the GWA is similar to HF in that the bare 
interaction in the exchange self-energy is replaced by a screened interaction and this screening 
is not sufficient to alleviate the deficiency of HF. This is an obvious problem when calculating 
molecular binding energies, and has been discussed in more detail in Ref. [21]. 

Let us now turn to calculations of atomic energy differences. It is evident from the shape 
of the binding curves around the equilibrium separation in Fig. |2 that SC-GW can produce 
accurate total energy differences. Calculations on atoms using the LW functional have also 
shown that two-electron removal energies, AE = En -2 — En, can be very accurately given 
within the GW approximation [11]. We therefore calculated the SC-GW removal energies of 
Be and Mg, as shown in Table [H] We find excellent agreement with the experimental results 
for both Mg and Be, the deviation being ten times smaller than those from the HF calculations. 
This improvement is in keeping with the results obtained by Shirley and Martin for GqWq 
calculations on atoms [4]. In Table ITTT1 we show the ionization potentials obtained from the 
EKT, both from the SC-GW and the non-selfconsistent GqWq Green function. The latter is 
obtained by iterating the Dyson equation once, starting from an LDA or HF Green function. 
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Table II - Two-electron removal energies En-2 — Em (in eV) calculated from SC-GW , compared to 
HF values and the experimental values. 



System 


SC-GW 


HF 


Expt. 1 


Mg - Mg 2+ 


22.59 


21.33 


22.68 


Be - Be 2+ 


27.59 


26.17 


27.53 



r From Ref. [22]. 



For most of the systems, the SC-GW ionization potentials are in good agreement with the 
experimental values, and in several cases better than those of GqWq. This is in contrast to 
the results for the electron gas, where self-consistency worsens the spectral properties [8]. 

The results for beryllium differ from those recently published by Delaney et al. [6]. We 
find a smaller difference between the SC-GW and the GoWo(LDA) results, and the latter 
value is also further away from the exact value than reported in Ref. [6]. One explanation 
for this deviation may be that while we obtained the ionization potentials from the EKT, 
Delaney et al. calculated them from the poles of the Fourier transformed function G(w). For 
the self-consistent ionization potentials, these methods should give the same result (they do 
in fact only differ with 0.2 eV), but for the GqWq Green function it is not obvious that the 
results should agree. Another difference is that we have carried out our calculations in a 
basis of Slater functions, while the orbitals in Ref. [6] are represented on a grid. The Slater 
basis was systematically extended until reaching convergence with respect to the total energy. 
We include HF orbitals with very large eigenenergies, e.g., for Be states up to 843 Hartree, 
while for Ne the highest orbital energy was 976 Hartree. We found good agreement between 
second-order M0ller-Plesset calculations with our basis sets and highly converged results from 
the literature [24]. This does not imply simultaneous convergence of other properties such as 
the ionization potential. In Table Hvl we illustrate the convergence of the beryllium atom for 
two different basis sets. The main difference between the sets is that basis I contains Slater 
functions optimized for HF calculations [23]. The uncertainty of ~ 0.02 eV in the ionization 
potential indicated in Table llVl is typical for the calculations on atoms presented in Table ITTT1 

Conclusions. - In summary, we have solved the Dyson equation within GWA to self- 
consistency for a number of atoms and diatomic molecules. We have shown that SC-GW 
gives good total energy differences and ionization potentials, significantly improving the HF 
results. We demonstrated that self-consistency improves the Go Wo ionization potentials for 
most systems studied and has the additional advantage of providing unambiguous results. 



Table III - Ionization potentials (eV) calculated from the EKT, using the self-consistent Green 
function and the Green function calculated from one iteration of the Dyson equation, starting from 
Glda and Ghf- 



System 


Go Wo (LDA) 


G Wo (HF) 


GW 


Expt. 1 


He 


23.65 


24.75 


24.56 


24.59 


Be 


8.88 2 


9.19 


8.66 2 


9.32 


Ne 


21.06 


21.91 


21.77 


21.56 


Mg 


7.52 


7.69 


7.28 


7.65 


H 2 


15.92 


16.52 


16.22 


15.43 


LiH 


6.87 


8.19 


7.85 


7.9 



^rom Ref. [22] 

2 To be compared with the Go Wo value 9.25 and the SC-GW value 8.47, reported in in Ref. [6] 
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Table IV - Convergence of the beryllium ionization potential, IP, (in eV) and total energy (in 
Hartrees) for two different basis sets. The value of l max indicates the maximum angular momentum 
quantum number used in the basis. 





/ — 9 

"max — ^ 


^max — 3 


^max — 4 


^max — 5 


^max — 6 


^max — 7 


IP: Basis I 


8.552 


8.602 


8.625 


8.636 


8.641 


8.644 


IP: Basis II 


8.439 


8.615 


8.637 


8.649 


8.654 


8.656 


E: Basis I 


-14.6954 


-14.6999 


-14.7016 


-14.7024 


-14.7028 


-14.7028 


E: Basis II 


-14.6807 


-14.6998 


-14.7015 


-14.7024 


-14.7027 


-14.7028 



Moreover, we have shown that the LW functional gives total energies in excellent agreement 
with the SC-GW energies, at a fraction of the computing time. This demonstrates the con- 
siderable usefulness of the LW functional for estimating the accuracy of various self-energy 
approximations . 

* * * 
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